7  Archaeobotany

Page under construction

The results presented here are preliminary and the chapter has yet to be written.

In this chapter, I will present the macrobotanical data from 190 case studies used to carry on this research (Chapter 3), along with the statistics performed on the data. The results will be first presented temporally, and a discussion of the diachronic trends will follow at the end of the chapter.

7.1 Case studies

The following map shows the sites under investigation, divided by chronology. Please select the desired chronology (or chronologies) from the legend at the top-right of the map.

Figure 7.1: Legend: R = Roman, LR = Late Roman, EMA = Early Middle Ages, Ma = 11th c. onwards

7.2 Ubiquity

In Chapter 6, ubiquity has been described as the best solution to present the archaeobotanical remains from the Italian peninsula, given the numerous biases in the samples. The heatmap below (Figure 7.2) provides a good overview of the temporal trends of presence of cereals, legumes, fruits and nuts in the entire area under examination.

Show the code

# Load the libraries
# Note: these libraries are used for the data visualizations in this page.
library(RColorBrewer)
library(reshape2)
library(ggplot2)
library(hrbrthemes)
library(plotly)
library(patchwork)

## UBIQUITY

## Creating a dataframe that contains the ubiquity of each century under examination. 
Ubiquity_table <- data.frame(
  "I BCE" = archaeobotany_tables(plants_export, -1)$Ubiquity_exp,  
  "I CE" = archaeobotany_tables(plants_export, 1)$Ubiquity_exp,
  "II CE" = archaeobotany_tables(plants_export, 2)$Ubiquity_exp,
  "III CE" = archaeobotany_tables(plants_export, 3)$Ubiquity_exp,
  "IV CE" = archaeobotany_tables(plants_export, 4)$Ubiquity_exp,
  "V CE" = archaeobotany_tables(plants_export, 5)$Ubiquity_exp,
  "VI CE" = archaeobotany_tables(plants_export, 6)$Ubiquity_exp,
  "VII CE" = archaeobotany_tables(plants_export, 7)$Ubiquity_exp,
  "VIII CE" = archaeobotany_tables(plants_export, 8)$Ubiquity_exp,
  "IX CE" = archaeobotany_tables(plants_export, 9)$Ubiquity_exp,
  "X CE" = archaeobotany_tables(plants_export, 10)$Ubiquity_exp,
  "XI CE" = archaeobotany_tables(plants_export, 11)$Ubiquity_exp
  )

# Transform the ubiquity table into a matrix
Ubiquity_mat <- as.matrix(Ubiquity_table) 

# Rename the centuries
colnames(Ubiquity_mat) <- c("1st c. BCE", "1st c. CE", "2nd c. CE",
                            "3rd c. CE", "4th c. CE", "5th c. CE",
                            "6th c. CE", "7th c. CE", "8th c. CE",
                            "9th c. CE", "10th c. CE", "11th c. CE") 

# The data has to be molten to use it with ggplot2
# (package: reshape2)
Ubiquity_melt <- melt(Ubiquity_mat)

# Let's now rename the columns 
colnames(Ubiquity_melt) <- c("Taxon", "Century", "Ubiquity")

# Add a column for the text tooltip
Ubiquity_melt <- Ubiquity_melt %>%
  mutate(text = paste0("Taxon: ", Taxon, "\n", "Century: ", Century, "\n", "Value: ",round(Ubiquity,2)))

# Create the heatmap with ggplot2
Ubiquity_HM <- ggplot(Ubiquity_melt, aes(Century, Taxon, fill=Ubiquity, text=text)) + 
  geom_tile(colour="white") +
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "right",
        axis.ticks = element_blank(), 
        axis.text.x = element_text(angle = 90, hjust = 0)
        ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Ubiquity",
    subtitle="Diachronical heatmap of recorded plant species"
  ) +
  scale_fill_gradient(low = "white", high = "black")

Figure 7.2: Diachronical heatmap of recorded plant species

7.2.1 Macroregional differences

The heatmap (Figure 7.2) shows the diachronical ubiquity values of the Italian mainland. It is possible however to compare plants’ ubiquity values in different areas of the peninsula. The R function Ubiquity_macroreg_chrono() (Section 2.3) was created to subset data from Northern, Central and Southern Italian regions, using modern boundaries. Archaeobotanical data from Italy is scarce, and subsetting the dataset required a larger chronological division to have enough sites for a valid statistical interpretation of the results. For this reason, the ubiquity values are presented using the variable Chronology rather than subsetting the individual centuries. For a clearer reading of the plot, the taxa have been divided into–Cereals, Pulses and Fruits/Nuts. Some taxa have been omitted from the plot.

Source code:

Show the code: data preparation

# Ubiquity by Italian Macro regions: Northern, Central and Southern Italy

# Load the libraries
library(vegan)
library(matrixStats)
library(patchwork)

# Creating a dataframe with the ubiquities of all macroregions and chronologies
bot_macroreg <- rbind(
  Ubiquity_R_NI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Northern Italy", "R"),
  Ubiquity_R_CI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Central Italy", "R"),
  Ubiquity_R_SI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Southern Italy", "R"),
  Ubiquity_LR_NI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Northern Italy", "LR"),
  Ubiquity_LR_CI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Central Italy", "LR"),
  Ubiquity_LR_SI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Southern Italy", "LR"),
  Ubiquity_EMA_NI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Northern Italy", "EMA"),
  Ubiquity_EMA_CI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Central Italy", "EMA"),
  Ubiquity_EMA_SI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Southern Italy", "EMA"),
  Ubiquity_Ma_NI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Northern Italy", "Ma"),
  Ubiquity_Ma_CI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Central Italy", "Ma"),
  Ubiquity_Ma_SI <- Ubiquity_macroreg_chrono(Archaeobot_Condensed,"Southern Italy", "Ma")
)

# Re-arranging the cereals/macroregions for visualisation on the Y axis
level_macroreg_order <- c("Southern Italy", "Central Italy", "Northern Italy")

level_cereals_order <- c("Common.Wheat", "Barley", "Rye", 
                         "Einkorn", "Emmer", "Proso.millet", 
                         "Foxtail.millet", "Oats", "Sorghum")

# Cereals
cer_ubiquity_macroreg.R <- filter(bot_macroreg, Chronology=="R" & Plant.Type=="Cereals")
cer_ubiquity_macroreg.R <- filter(cer_ubiquity_macroreg.R, Macroregion!="Central Italy")
cer_ubiquity_macroreg.LR <- filter(bot_macroreg, Chronology=="LR" & Plant.Type=="Cereals")
cer_ubiquity_macroreg.EMA <- filter(bot_macroreg, Chronology=="EMA" & Plant.Type=="Cereals")
cer_ubiquity_macroreg.Ma <- filter(bot_macroreg, (Chronology=="Ma" & Plant.Type=="Cereals"))
cer_ubiquity_macroreg.Ma <- filter(cer_ubiquity_macroreg.Ma, Macroregion!="Southern Italy")

#Pulses
puls_ubiquity_macroreg.R <- filter(bot_macroreg, Chronology=="R" & Plant.Type=="Pulses")
puls_ubiquity_macroreg.R <- filter(puls_ubiquity_macroreg.R, Macroregion!="Central Italy")
puls_ubiquity_macroreg.R <- filter(puls_ubiquity_macroreg.R, Plant!="Chickpea")
puls_ubiquity_macroreg.LR <- filter(bot_macroreg, Chronology=="LR" & Plant.Type=="Pulses")
puls_ubiquity_macroreg.LR <- filter(puls_ubiquity_macroreg.LR, 
                                    Macroregion!="Southern Italy")
puls_ubiquity_macroreg.LR <- filter(puls_ubiquity_macroreg.LR, Plant!="Chickpea")
puls_ubiquity_macroreg.EMA <- filter(bot_macroreg, Chronology=="EMA" & Plant.Type=="Pulses")
puls_ubiquity_macroreg.Ma <- filter(bot_macroreg, Chronology=="Ma"  & Plant.Type=="Pulses")
puls_ubiquity_macroreg.Ma <- filter(puls_ubiquity_macroreg.Ma, 
                                    Macroregion!="Southern Italy")

#Fruits (+ Subset)
fnuts_ubiquity_macroreg.R <- filter(bot_macroreg, Chronology=="R" & Plant.Type=="Fruits/Nuts")
fnuts_ubiquity_macroreg.R <- subset(fnuts_ubiquity_macroreg.R, (Plant == "Wild.Cherry" | Plant == "Walnut" | Plant == "Peach" | Plant == "Olive" |Plant == "Grape" | Plant =="Fig" | Plant =="Apple"))
fnuts_ubiquity_macroreg.R <- filter(fnuts_ubiquity_macroreg.R, Macroregion!="Central Italy")
fnuts_ubiquity_macroreg.LR <- filter(bot_macroreg, Chronology=="LR" & Plant.Type=="Fruits/Nuts")
fnuts_ubiquity_macroreg.LR <- subset(fnuts_ubiquity_macroreg.LR, (Plant == "Wild.Cherry" | Plant == "Walnut" | Plant == "Peach" | Plant == "Olive" |Plant == "Grape" | Plant =="Fig" | Plant =="Apple"))
fnuts_ubiquity_macroreg.EMA <- filter(bot_macroreg, Chronology=="EMA" & Plant.Type=="Fruits/Nuts")
fnuts_ubiquity_macroreg.EMA <- subset(fnuts_ubiquity_macroreg.EMA, (Plant == "Wild.Cherry" | Plant == "Walnut" | Plant == "Peach" | Plant == "Olive" |Plant == "Grape" | Plant =="Fig" | Plant =="Apple"))
fnuts_ubiquity_macroreg.Ma <- filter(bot_macroreg, Chronology=="Ma"  & Plant.Type=="Fruits/Nuts")
fnuts_ubiquity_macroreg.Ma <- subset(fnuts_ubiquity_macroreg.Ma, (Plant == "Wild.Cherry" | Plant == "Walnut" | Plant == "Peach" | Plant == "Olive" |Plant == "Grape" | Plant =="Fig" | Plant =="Apple"))
fnuts_ubiquity_macroreg.Ma <- filter(fnuts_ubiquity_macroreg.Ma, Macroregion!="Southern Italy")

Show the code: plots

# Cereals plots ubiquity
cer_ubiquity_macroreg_R.HM <- ggplot(cer_ubiquity_macroreg.R, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  factor(Plant, levels=rev(level_cereals_order)),
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Roman"
  ) +  scale_fill_gradient(low = "white", high = "black")


cer_ubiquity_macroreg_LR.HM <- ggplot(cer_ubiquity_macroreg.LR, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  factor(Plant, levels=rev(level_cereals_order)),
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Late Roman"
  ) +  scale_fill_gradient(low = "white", high = "black")

cer_ubiquity_macroreg_EMA.HM <- ggplot(cer_ubiquity_macroreg.EMA, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  factor(Plant, levels=rev(level_cereals_order)),
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Early Medieval"
  ) +  scale_fill_gradient(low = "white", high = "black")

cer_ubiquity_macroreg_Ma.HM <- ggplot(cer_ubiquity_macroreg.Ma, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  factor(Plant, levels=rev(level_cereals_order)),
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Medieval"
  ) +  scale_fill_gradient(low = "white", high = "black")


Cereals_Ubiquity_MacroReg_Patchwork <- (cer_ubiquity_macroreg_R.HM|cer_ubiquity_macroreg_LR.HM)/(cer_ubiquity_macroreg_EMA.HM|cer_ubiquity_macroreg_Ma.HM)
Cereals_Ubiquity_MacroReg_Patchwork + plot_annotation(
  title = 'Cereals',
  subtitle = 'Ubiquity (%), plotted by macroregion and chronology.',
  caption='Note: Data was too scarce for Roman Central Italy and Medieval Southern Italy.'
  )

Show the code: plots

# Pulses plots ubiquity
puls_ubiquity_macroreg_R.HM <- ggplot(puls_ubiquity_macroreg.R, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="#cfcfcf", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Roman"
  ) +  scale_fill_gradient(low = "white", high = "black")


puls_ubiquity_macroreg_LR.HM <- ggplot(puls_ubiquity_macroreg.LR, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="#ffffff", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Late Roman"
  ) +  scale_fill_gradient(low = "white", high = "black")

puls_ubiquity_macroreg_EMA.HM <- ggplot(puls_ubiquity_macroreg.EMA, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Early Medieval"
  ) +  scale_fill_gradient(low = "white", high = "black")

puls_ubiquity_macroreg_Ma.HM <- ggplot(puls_ubiquity_macroreg.Ma, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Medieval"
  ) +  scale_fill_gradient(low = "white", high = "black")


Pulses_Ubiquity_MacroReg_Patchwork <- (puls_ubiquity_macroreg_R.HM|puls_ubiquity_macroreg_LR.HM)/(puls_ubiquity_macroreg_EMA.HM|puls_ubiquity_macroreg_Ma.HM)
Pulses_Ubiquity_MacroReg_Patchwork + plot_annotation(
  title = 'Pulses',
  subtitle = 'Ubiquity (%), plotted by macroregion and chronology.',
  caption='Note: Data was too scarce for Roman Central Italy and Late Roman/Medieval Southern Italy.'
)

Show the code: plots

# Fruits nuts plots

fnuts_ubiquity_macroreg_R.HM <- ggplot(fnuts_ubiquity_macroreg.R, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Roman"
  ) +  scale_fill_gradient(low = "white", high = "black")


fnuts_ubiquity_macroreg_LR.HM <- ggplot(fnuts_ubiquity_macroreg.LR, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Late Roman"
  ) +  scale_fill_gradient(low = "white", high = "black")

fnuts_ubiquity_macroreg_EMA.HM <- ggplot(fnuts_ubiquity_macroreg.EMA, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Early Medieval"
  ) +  scale_fill_gradient(low = "white", high = "black")

fnuts_ubiquity_macroreg_Ma.HM <- ggplot(fnuts_ubiquity_macroreg.Ma, aes(
  factor(Macroregion, levels=(level_macroreg_order)),
  Plant,
  fill=(Ubiquity)
)) + 
  geom_tile(colour="white") +
  geom_text(aes(label = Ubiquity), colour="white", size=3)+ 
  scale_alpha(range=c(0,1)) +
  scale_x_discrete("", expand = c(0, 0)) + 
  scale_y_discrete("", expand = c(0, 0)) + 
  theme_grey(base_size = 9) + 
  theme(legend.position = "none",
        axis.ticks = element_blank()
  ) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+
  labs(
    title="Medieval"
  ) +  scale_fill_gradient(low = "white", high = "black")


FrNuts_Ubiquity_MacroReg_Patchwork <- (fnuts_ubiquity_macroreg_R.HM|fnuts_ubiquity_macroreg_LR.HM)/(fnuts_ubiquity_macroreg_EMA.HM|fnuts_ubiquity_macroreg_Ma.HM)
FrNuts_Ubiquity_MacroReg_Patchwork + plot_annotation(
  title = 'Fruits/Nuts',
  subtitle = 'Ubiquity (%), plotted by macroregion and chronology.',
  caption='Note: Data was too scarce for Roman Central Italy and Medieval Southern Italy.'
)

7.2.1.1 Cereals

It is interesting to notice how in the Roman age, cereals are similarly ubiquitous in Southern and Northern Italy, although there are some exceptions (i.e. einkorn, rye, oats, and proso millet) that can derive from the randomness of sampling. Unfortunately, only three sites provided botanical samples for Roman Central Italy and their values have then been omitted from the plot. These sites (all from Tuscany) were studied by the Roman Peasant Project (Bowes and University of Pennsylvania, 2020) and only reported three kinds of cereal: common wheat, emmer, and barley. Similar ubiquity values from Northern and Southern Italy in the Roman age may suggest similar production patterns in the whole Italian mainland, even though more data is required. In the Late Roman age, ubiquity data has been calculated for the three macroregions, with Southern Italy being the least trustworthy (five sites in total). Three crops are found on 62.5-75% of the Central Italian sites: common wheat, barley and emmer. Other cereals are present, but less ubiquitously. The cereal triad aforementioned seems to be diffused in the south as well. Conversely, in Northern Italy common wheat and barley were indeed important cultivations but had to compete with other cereals including millet, sorghum, and rye. The latter had a significant increase, being present on almost 30% of the Northern sites (as opposed to the Roman 15%). The Early Medieval age seems to mark a shift in Italian agricultural practices, as cereal ubiquities are much more variable regionally. In Southern Italy, the triad of common wheat, barley and emmer were still the predominant cereals. These cereals are ubiquitous in Central and Northern Italy as well, although these regions adopt polyculture with a diversified number of cereals. The samples from the Medieval age are fewer in number since the upper boundary of this project’s chronology is the 11^th c. Despite the short chronology, it is possible to make some considerations. Medieval Centraly Italy relied heavily on common wheat, barley and emmer, with other cereals increasingly important. Barley is the most ubiquitous cereal in Northern Italy in this period, followed by common wheat, millets and sorghum.

Figure 7.3: Diachronical heatmap of cereals in the Italian macroregions

7.2.1.2 Pulses

In the Roman Age, pulses are an important part of the diet and are cultivated both in Northern and Southern Italy. In the latter, vetch/broad beans are present in 22-32% of the samples, and lentils are present in 38% of the sites. In the Late Roman Age, broad beans are equally important in Central and Northern Italy, and peas are present in 50% of the Central Italian sites. In the Early Medieval Age, pulses are present in many Central Italian sites, especially blue/red peas, broad beans and other Fabaceae. Lentils and broad beans are also cultivated in almost half of the Northern Italian sites. The importance of pulses in Central Italy is confirmed by the 11^th c. samples, where every specie is present in over 66% of the sites and Fabaceae and blue/red peas are found in every sample. Conversely, in Northern Italy broad bean is found in 66% of the sites.

Figure 7.4: Diachronical heatmap of pulses in the Italian macroregions

7.2.1.3 Fruits and nuts

Olive and grape are two essential cultivations in the Italian peninsula. Olive pits, as can be expected, are more ubiquitous in Southern Italy, where in Roman times are present in >87% of the sites and in over 58% of the sites in the following chronologies¹. Conversely, grape is important in Central and Northern Italy in the Late Roman, Early Medieval and Medieval ages.

Figure 7.5: Diachronical heatmap of fruits/nuts in the Italian macroregions

7.3 Richness and diversity

Section in progress

7.3.1 Richness and diversity in the Italian macroregions

Show the code: data preparation

# Species richness based on geographical features
# RELATIVE PROPORTIONS OF ARCHAEOBOT_VIZ QUERY EXPORT FROM THE DB
# (Condensed, without totals)

# Remove NAs
Df_Cond_Plants[is.na(Df_Cond_Plants)] <-0 

# Generate a dataframe with the relative proportions and round the results
Df_Cond_Plants_Rel <- decostand(Df_Cond_Plants[11:50], method = "total")
Df_Cond_Plants_Rel <- round(Df_Cond_Plants_Rel, digits=2)

# Add more info to the dataframe
Df_Cond_Plants_Rel_Richness_Diversity <- data.frame(
                                 "Geo" = Df_Cond_Plants$Geo,
                                 "Chronology" = Df_Cond_Plants$Chronology, 
                                 "Type"= Df_Cond_Plants$Type, 
                                 "Specnumber" = specnumber(Df_Cond_Plants_Rel),
                                 "Shannon Div" = diversity(Df_Cond_Plants_Rel),
                                 Df_Cond_Plants_Rel
)

Df_Cond_Plants_Rel_withMacroregion <- data.frame("Geo" = Df_Cond_Plants$Geo,
                                                 "Chronology" = Df_Cond_Plants$Chronology,
                                                 "Type"= Df_Cond_Plants$Type, 
                                                 "Macroregion" = Df_Cond_Plants$name_macroreg,
                                                 "Specnumber" = specnumber(Df_Cond_Plants_Rel[1:10]), #Only cereals
                                                 "Shannon Div" = diversity(Df_Cond_Plants_Rel[1:10]),
                                                 Df_Cond_Plants_Rel[1:10]
)

# Let's plot the diversity by macroregion

# Creating the dataframes for R and EMA age 
# I know it's called "Plants" but it's actually just cereals
Df_Cond_Plants_Rel_withMacroregion.R <- filter(Df_Cond_Plants_Rel_withMacroregion, Chronology == "R")
Df_Cond_Plants_Rel_withMacroregion.LR <- filter(Df_Cond_Plants_Rel_withMacroregion, Chronology == "LR")
Df_Cond_Plants_Rel_withMacroregion.EMA <- filter(Df_Cond_Plants_Rel_withMacroregion, Chronology == "EMA")

Show the code: plots

pal_RichnessvsGeo <- c("cadetblue3", "gold1",  "bisque4", "palegreen4")

plot_RichnessMacroReg.R <- ggplot(Df_Cond_Plants_Rel_withMacroregion.R, aes(x = Macroregion, y = Specnumber, fill = Macroregion)) +
  geom_violin(trim=FALSE) + 
  geom_boxplot(width=0.1, fill="white")+
  scale_fill_manual(values = pal_RichnessvsGeo) +
  geom_jitter(alpha=0.3)+
  scale_x_discrete(labels = c("Central Italy \n (n = 5)", "Northern Italy \n (n = 41)", "Southern Italy \n (n=34)")) +
  theme(legend.position = "none",
        plot.background = element_rect("white"),
        panel.background = element_rect("white"),
        panel.grid = element_line("grey90"),
        axis.line = element_line("gray25"),
        axis.text = element_text(size = 12, color = "gray25"),
        axis.title = element_text(color = "gray25"),
        legend.text = element_text(size = 12)) + 
  labs(x = "Macroregion",
       y = "Number of species per site",
       title = "R - Cereal richness")

plot_RichnessMacroReg.EMA <- ggplot(Df_Cond_Plants_Rel_withMacroregion.EMA, aes(x = Macroregion, y = Specnumber, fill = Macroregion)) +
  geom_violin(trim=FALSE) + 
  geom_boxplot(width=0.1, fill="white")+
  scale_fill_manual(values = pal_RichnessvsGeo) +
  geom_jitter(alpha=0.3)+
  scale_x_discrete(labels = c("Central Italy \n (n = 11)", "Northern Italy \n (n = 36)", "Southern Italy \n (n=17)")) +
  theme(legend.position = "none",
        plot.background = element_rect("white"),
        panel.background = element_rect("white"),
        panel.grid = element_line("grey90"),
        axis.line = element_line("gray25"),
        axis.text = element_text(size = 12, color = "gray25"),
        axis.title = element_text(color = "gray25"),
        legend.text = element_text(size = 12)) + 
  labs(x = "Macroregion",
       y = "Number of species per site",
       title = "EMA - Cereal richness")

# Check median of cereal richness per macroregion

# ROMAN NORTH
mean((filter(Df_Cond_Plants_Rel_withMacroregion.R, Macroregion=="Northern Italy"))$Specnumber)
#ROMAN SOUTH
mean((filter(Df_Cond_Plants_Rel_withMacroregion.R, Macroregion=="Southern Italy"))$Specnumber)
# EMA NORTH
mean((filter(Df_Cond_Plants_Rel_withMacroregion.EMA, Macroregion=="Northern Italy"))$Specnumber)
#EMA SOUTH
mean((filter(Df_Cond_Plants_Rel_withMacroregion.EMA, Macroregion=="Southern Italy"))$Specnumber)

Cereals share similar presence values in Roman Northern and Southern Italian sites (Figure 7.6). Central Italy reports higher values, although this is based only on three sites and hence it is not reliable. During the Early Middle Ages, Central Italy again is the richest in cereals, closely followed by Northern Italy. Interestingly, Southern Italy still reports values very close to the Roman age. A full list of the Southern Italian EMA sites is reported in Table 7.1.

(a) Roman age.

(b) Early Medieval age.

Figure 7.6: Violin plots of cereal richness in the Italian macroregions. The grey dots (jitters) indicate the value for the single site, while the white boxplot shows the median and the quartile values.

Table 7.1: List of Southern Italian sites with chronology EMA

ID	Site	Region	Geography	Type	Culture/Influence
98	S. Maria in Cività, D85	Molise	Hilltop	Urban	Lombard
107	S. Giovanni di Ruoti, Phase 3A	Basilicata	Mountain	Monastery	Lombard
107	S. Giovanni di Ruoti, Phase 3B	Basilicata	Mountain	Monastery	Lombard
198	Salapia, area botteghe, US 2475	Puglia	Coast/Lagoon	Urban	Lombard
198	Salapia, area botteghe, US 2437	Puglia	Coast/Lagoon	Urban	Lombard
199	Salapia, area conceria, US 2054	Puglia	Coast/Lagoon	Urban	Lombard
199	Salapia, area conceria, US 2211-2217	Puglia	Coast/Lagoon	Urban	Lombard
199	Salapia, area conceria, 8th-9th c.	Puglia	Coast/Lagoon	Urban	Lombard
196	Faragola, wastepit 61	Puglia	Plain	Rural, villa	Lombard
196	Faragola, wastepit 66	Puglia	Plain	Rural, villa	Lombard
234	Colle Castellano, Phase 3-4	Molise	Hill	Urban	Lombard
177	San Vincenzo al Volturno, kitchen area	Molise	Hill	Monastery	Lombard
101	Supersano, loc. Scorpo	Puglia	Plain	Rural	Byzantine
250	Apigliano, 9th-10th c., pits	Puglia	Plain	Rural	Byzantine
250	Apigliano, 10th-11th c., pits	Puglia	Plain	Rural	Byzantine
196	Faragola, granary A7	Puglia	Plain	Rural, villa	Lombard
196	Faragola, granary A8	Puglia	Plain	Rural, villa	Lombard

7.4 Cereals regionality

7.4.1 PERMANOVA

Notes on terminology: PERMANOVA

Permutational multivariate analysis of variance (PERMANOVA) is a non-parametric multivariate statistical test used to compare group of objects. By using measure space, the null hypothesis that the centroids and dispersion of groups are identical is tested. The null hypothesis is rejected if either the centroid or the spread of the objects differs between the groups. A prior calculation of the distance between any two objects included in the experiment is used to determine whether the test is valid or not² (Anderson, 2017). In this context, the null hypothesis is that there is no regional difference in the cereals dataset, with cereals being evenly distributed across macroregions and chronologies.

The suggestion of an Early Medieval shift in cereal farming stated in Section 7.2.1 and Section 7.3.1 needs statistical support. Considering that data is not unimodal and that we are dealing with presence/absence analysis, the best choice is to use a non-parametric test as PERMANOVA on the early medieval botanical dataset. Prior to performing the test, it was necessary to pre-process data by:

• Selecting all the cereals columns of the plant remains table, keeping the categorical variables: Macroregion.

• Removing empty rows (samples that only reported seeds/fruits, but not cereals).

• Transforming the raw counts into presence/absence, using the function decostand() (method=pa) in the R package vegan (Oksanen et al., 2020).

Show the code: Pre-processing

# Testing the results: Regionality in the dataset? 
library(vegan)

set.seed(29)

# Pre processing: remove empty rows

# Note: The input table is the CONDENSED table without totals

# Selecting all the cereals columns of the plant remains table, keeping some categorical variables 
cer_macroreg_ubiquity_transp.tot <- Df_Cond_Plants[c(4,5,6,7,11:19)]

# Selecting all rows with data (since we selected only with cereals, and some sites only had fruits/pulses we might have empty rows)
cer_macroreg_ubiquity_transp.tot <- cer_macroreg_ubiquity_transp.tot[rowSums(cer_macroreg_ubiquity_transp.tot[5:13])>0,]

# Assigning a column name "Macroregion"
colnames(cer_macroreg_ubiquity_transp.tot)[1] = "Macroregion"

# Selecting the Chronology of interest (EMA) and excluding Central Italy
cer_macroreg_ubiquity_transp.tot<- filter(cer_macroreg_ubiquity_transp.tot, Macroregion!="Central Italy" & Chronology=="EMA")

# dividing categorical and numerical columns
cer_macroreg_ubiquity_transp.categ <- cer_macroreg_ubiquity_transp.tot[1:4]
cer_macroreg_ubiquity_transp.data <- cer_macroreg_ubiquity_transp.tot[5:13]

# Converting the numerical columns into a presence/absence matrix (using method=pa)
cer_macroreg_ubiquity_transp.dist <- decostand(cer_macroreg_ubiquity_transp.data, method="pa", na.rm=TRUE)

After the pre-processing, it was possible to run the PERMANOVA using the function adonis2() in the package vegan. The function creates a distance matrix and computes an analysis of variance on the matrix. The method chosen to calculate the distance matrix is the jaccard distance. The Jaccard distance, based on the Jaccard similarity index, is a value of dissimilarity between sample sets (Kosub, 2019). When compared to other dissimilarity indices, it is more appropriate for presence/absence analyses as it is not based on Euclidean distance.

Results of adonis2().

Show the code: adonis2()

cer_macroreg_ubiquity_transp.div <- adonis2(
  cer_macroreg_ubiquity_transp.dist ~ Macroregion, 
  data = cer_macroreg_ubiquity_transp.categ,
  permutations = 10000, method="jaccard"
  )

Permutation test for adonis under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 10000

adonis2(formula = cer_macroreg_ubiquity_transp.dist ~ Macroregion, data = cer_macroreg_ubiquity_transp.categ, permutations = 10000, method = "jaccard")
            Df SumOfSqs      R2      F    Pr(>F)    
Macroregion  1   1.3554 0.13066 7.8153 9.999e-05 ***
Residual    52   9.0181 0.86934                     
Total       53  10.3734 1.00000                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The calculation of PERMANOVA on the Roman dataset using the variable Macroregion reported a not statistically significant p-value, suggesting a degree of homogeneity. On the other hand, when applied to the early Medieval dataset, the results were highly significant (0<p<0.001) with a 99.99% confidence. PERMANOVA has not been calculated on the Late Roman dataset as samples from Southern Italy were too scarce. To confirm the significance of the results of PERMANOVA applied to the EMA dataset, it is necessary to check that its assumptions are met (especially since we are dealing with small groups of data). Firstly, we check the homogeneity of variances. The function betadisper() from the package vegan provides the distances of group samples from centroids. If the variation is even, the null hypothesis of no difference in dispersion between groups is accepted. To test the variation, it is possible to use the analysis of variance (ANOVA).

Show the code: betadisper()

# We do not need to calculate the distance separately, but it will be useful later for the betadisper() function

# Distance dissimilarity matrix with the Jaccard method
cer_macroreg_ubiquity_transp.dist2 <- vegdist(cer_macroreg_ubiquity_transp.dist, method="jaccard", na.rm=TRUE)

# Betadisper: distances of group samples from centroids
cer_macroreg_ubiquity_transp.betadisper <- betadisper(cer_macroreg_ubiquity_transp.dist2, cer_macroreg_ubiquity_transp.categ$Macroregion)

Results of anova() on the betadisper.

Show the code: ANOVA on betadisper()

# We will see that the ANOVA's p-value is not significant meaning that group dispersions are homogenous 
#("Null hypothesis of no difference in dispersion between groups"; https://www.rdocumentation.org/packages/vegan/versions/2.4-2/topics/betadisper).

anova(cer_macroreg_ubiquity_transp.betadisper) # This should not be significant!

Analysis of Variance Table

Response: Distances
          Df Sum Sq   Mean Sq F value Pr(>F)
Groups     1 0.0055 0.0054954  0.1752 0.6773
Residuals 52 1.6311 0.0313675               

(a) Groups dispersions plot with confidence ellipses.

(b) Boxplot showing equal distances from centroid.

Figure 7.7: Results of the betadisper() (groups dispersions) on the distance matrix calculated with the Jaccard method.

The betadisper() graphs (Figure 7.7) show similar distances from the centroids for the categories Northern Italy and Southern Italy. In addition, the ANOVA on the betadisper() shows that the separation is not significant (p-value over the significance threshold), meaning that the groups dispersions are homogeneous. We can thus be confident in the PERMANOVA results and accept the difference between the two groups of sites under investigation. In other words, the Southern and Northern Italian group of sites are different during the Early Middle Ages for what concerns cereal farming.

Running the same tests on the Roman sites failed to separate the two groups of sites, confirming that there was not a significant difference in the types of cereals cultivated during the Roman age between Northern and Southern Italy.

7.4.2 nMDS

Notes on terminology: Wasserstein metric

The Wasserstein distance (or earth’s mover distance) is a measure of distance between two probability distributions on a metric space.

In addition to statistically testing the separation between the Northern and Southern Italian early medieval cereals dataset (Section 7.4.1), it is possible to measure the distance between groups of sites both in the Roman and early Middle ages. For this task, a dimensionality reduction algorithm has been chosen: the non-metric multidimensional scaling, from the R library vegan. A more in-depth explanation of the algorithm can be read in Section 6.6.2.3. To avoid fallacy in computations, the macroregion Central Italy and the chronologies LR (Late Roman) and Ma (11^th c. onwards) have been excluded from this test—the uneven distribution of the group of samples required a cautious approach. The nMDS has been run with a reduction to only one dimension, using KDE plots to visualize the results. Setting the dimension to one allows easier calculations of distance. In Figure 7.8 (a), it is possible to see the nMDS performed on the Roman cereals presence/absence dataset. As already pointed out, the PERMANOVA did not produce significant results for this dataset and the Wasserstein distance (calculated with the wasserstein1d() function in the transport library) is indeed shorter for the Roman dataset. For both chronologies there is an overlap in the curves, which is more considerable in the Roman age (indicating that the group of samples are more similar). The overlap for the EMA groups (Figure 7.8, b) is probably due to the fact that the presence of the noble grains is not by itself a ‘marker’ of Southern Italian sites—these grains are also very common in the North. The difference is that in the South noble grains are not cultivated in conjunction with other grains. The graph for the EMA chronology shows a clearer separation of the macroregional groups, with some minor overlaps. Moreover, the graph also displays variability in the Northern Italian dataset. The variability can also be assessed from the outliers in the boxplots in Figure 7.9.

Source code:

Show the code: Data preparation

# DATA PREPARATION #

# Creating a simple dataframe just for the purpose of the NCA
Df_Cond_Plants_simplified <- data.frame(
  "Chronology" = Df_Cond_Plants$Chronology, 
  "Type" = Df_Cond_Plants$Type,
  "Macroregion" = Df_Cond_Plants$name_macroreg,
  "Weight" = Df_Cond_Plants$weight,
  Df_Cond_Plants[11:ncol(Df_Cond_Plants)]
  )

# Joining site typology categories to simplify the output, using library stringr
library(vegan)
library(stringr)

Df_Cond_Plants_simplified$Type <- str_replace_all(
  Df_Cond_Plants_simplified$Type, c(
    "Rural site, villa" = "Rural", 
    "Rural site, mansio" = "Rural",
    "Religious, monastery" = "Religious",
    "Castle" = "Fortified",
    "Castrum" = "Fortified")
  )

# Convert to Presence (1) / Absence (0) using library vegan
Df_Cond_Plants_simplified[5:ncol(Df_Cond_Plants_simplified)] <- decostand(Df_Cond_Plants_simplified[5:ncol(Df_Cond_Plants_simplified)], method="pa")

# Creating dataframes for the chronologies R and EMA
Df_Cereals_R <- filter(Df_Cond_Plants_simplified, Chronology=="R" & Macroregion!="Central Italy")
Df_Cereals_EMA <- filter(Df_Cond_Plants_simplified, Chronology=="EMA" & Macroregion!="Central Italy")

# Selecting only cereals from both dataframes
Df_Cereals_R <- Df_Cereals_R[1:14]
Df_Cereals_R[is.na(Df_Cereals_R)] <- 0
Df_Cereals_R <- Df_Cereals_R[rowSums(Df_Cereals_R[,5:14] != 0) > 0, ] # Removing empty rows

Df_Cereals_EMA <- Df_Cereals_EMA[1:14]
Df_Cereals_EMA[is.na(Df_Cereals_EMA)] <- 0
Df_Cereals_EMA <- Df_Cereals_EMA[rowSums(Df_Cereals_EMA[,5:14] != 0) > 0, ]

(a) Roman age.

(b) Early Medieval Age.

Figure 7.8: nMDS reduction to one dimension performed on the Roman and early Medieval cereal datasets. The dashed colored lines indicate the mean of the distribution.

(a) Roman age.

(b) Early Medieval age.

Figure 7.9: Boxplots showing the nMDS scores for Northern and Southern Roman Italy.

Old Python code for the NCA to delete:

Show the code: Libraries

# Load Python libraries
#!pip install pandas
import pandas as pd
import os

#!pip install scikit-hubness
import random
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from sklearn.neighbors import (KNeighborsClassifier,
                                 NeighborhoodComponentsAnalysis)

from scipy.stats import wasserstein_distance

import seaborn as sns

# Set seed
random.seed(10)

Show the code: Selecting random samples, with replacement

#R
df_R_SI = df_R[df_R["Macroregion"]=="Southern Italy"].sample(20, random_state=21, replace="TRUE")
df_R_NI = df_R[df_R["Macroregion"]=="Northern Italy"].sample(20, random_state=21, replace="TRUE")
# Create a dataset with northern and southern Italy
df_R_merge = pd.concat([df_R_SI, df_R_NI], ignore_index=True)

data_R_byPlantGroup = df_R_merge.drop(['Chronology','Type', 'Macroregion', 'Weight'], axis=1)
labels_R_byPlantGroup = df_R_merge.iloc[:,2] # nrow, 0 for Chronology - nrow, 1 for Type - nrow,2 for Macroregion

#EMA
df_EMA_SI = df_EMA[df_EMA["Macroregion"]=="Southern Italy"].sample(20, random_state=21, replace="TRUE")
df_EMA_NI = df_EMA[df_EMA["Macroregion"]=="Northern Italy"].sample(20, random_state=21, replace="TRUE")
df_EMA_merge = pd.concat([df_EMA_SI, df_EMA_NI], ignore_index=True)
 
data_EMA_byPlantGroup = df_EMA_merge.drop(['Chronology','Type', 'Macroregion', 'Weight'], axis=1)
labels_EMA_byPlantGroup = df_EMA_merge.iloc[:,2] # nrow, 0 for Chronology - nrow, 1 for Type - nrow,2 for Macroregion

Show the code: Performing the NCA

#R
nca_for_KDE_R_PlantGroup = NeighborhoodComponentsAnalysis(n_components =1, init="lda").fit(data_R_byPlantGroup, labels_R_byPlantGroup)
reduction_for_KDE_R_PlantGroup = nca_for_KDE_R_PlantGroup.transform(data_R_byPlantGroup)
df_R_merge["value"] = reduction_for_KDE_R_PlantGroup
     
#EMA
nca_for_KDE_EMA_PlantGroup = NeighborhoodComponentsAnalysis(n_components =1, init="lda").fit(data_EMA_byPlantGroup, labels_EMA_byPlantGroup)
reduction_for_KDE_EMA_PlantGroup = nca_for_KDE_EMA_PlantGroup.transform(data_EMA_byPlantGroup)
df_EMA_merge["value"] = reduction_for_KDE_EMA_PlantGroup

Show the code: Wasserstein distance

# R
df_R_North = df_R_merge[df_R_merge["Macroregion"] == "Northern Italy"]
df_R_South = df_R_merge[df_R_merge["Macroregion"] == "Southern Italy"]

print('Roman age - Wasserstein Distance')
wasserstein_distance(df_R_South["value"], df_R_North["value"], u_weights=df_R_North["Weight"], v_weights=df_R_South["Weight"])


#EMA
print('EMA - Wasserstein Distance')
df_EMA_North = df_EMA_merge[df_EMA_merge["Macroregion"] == "Northern Italy"]
df_EMA_South = df_EMA_merge[df_EMA_merge["Macroregion"] == "Southern Italy"]

wasserstein_distance(df_EMA_South["value"], df_EMA_North["value"], u_weights=df_EMA_North["Weight"], v_weights=df_EMA_South["Weight"])

Show the code: Plotting the NCA

NCA_KDE_1D, ax = plt.subplots(1, 2, figsize=(10, 5), sharey=True, sharex=True)

sns.kdeplot(data=df_R_merge, x="value", ax=ax[0], hue="Macroregion", fill=True, alpha=.1, palette="colorblind", linewidth=1, legend=None).set(xlabel='NCA', title="(a) Roman age")

sns.kdeplot(data=df_EMA_merge, x="value", ax=ax[1], hue="Macroregion", fill=True, alpha=.1, palette="colorblind", linewidth=1).set(xlabel='NCA', title="(b) early Middle ages")

NCA_KDE_1D.text(0.08, 0.03, 'Weighted Wasserstein Distance: W = 20.67 \nPERMANOVA test: p>0.05', fontsize=10)
NCA_KDE_1D.text(0.68, 0.03, 'Weighted Wasserstein Distance: W = 88.32\nPERMANOVA test: p<0.001', fontsize=10)
plt.tight_layout()
plt.subplots_adjust(bottom=0.19)
plt.show()

7.4.3 GLM: Minor grains

Mapping the use of minor grains by chronology in each of the three macroregions.

Code

Minor_Grains <- with(
  Df_Cond_Plants.pa,
  data.frame(
  Chronology = Chronology,
  From.Century=From.Century, 
  To.Century=To.Century,
  Macroregion = name_macroreg,
  Site_Type = as.factor(Type),
  Geo = Geo,
  Grain_Count = Einkorn+Oats+Rye+Proso.millet+Foxtail.millet+Sorghum+Emmer,
  All_Grain_Count = rowSums(Df_Cond_Plants.pa[,c(11:19)])
  )
)

#Simplify categories
Minor_Grains$Site_Type <- str_replace(Minor_Grains$Site_Type, "Religious, monastery", "Religious")    
Minor_Grains$Site_Type <- str_replace(Minor_Grains$Site_Type, "Castle", "Fortified")    
Minor_Grains$Site_Type <- str_replace(Minor_Grains$Site_Type, "Castrum", "Fortified")  

# Transform categories to factors
Minor_Grains[,c(1,4,5,6)] <- lapply(Minor_Grains[,c(1,4,5,6)], factor)

# Count the occurrences of minor grains in each context
Minor_Grains$Present <- ifelse(Minor_Grains$Grain_Count > 0, 1, 0)

# Test: Transform Chronology to discrete values so to visualise 
Minor_Grains <- Minor_Grains %>% 
    mutate(Century = map2(From.Century, To.Century, `:`)) %>% 
    select(-From.Century, -To.Century) %>% 
    unnest
# End

# Remove duplicate rows caused by uncertain chronology
#Minor_Grains <- Minor_Grains[!duplicated(Minor_Grains[-1]), ]

# Create a copy of the dataframe
Minor_Grains2 <- Minor_Grains

# Convert to list
Minor_Grains2.list <- list(
    MG_Found = Minor_Grains2$Grain_Count,
    MG_Tot = 7,
    RC_ID = as.factor((interaction(Minor_Grains2$Chronology, 
                                   Minor_Grains2$Macroregion)))
   )

# Create a binomial model, where N is 7 because that is the total number of minor 
# grains studied in this thesis
m_mg_reg <- ulam(
    alist(
        MG_Found ~ dbinom( 7 , p ),
        logit(p) <- RegChr[RC_ID],
        RegChr[RC_ID] ~ dnorm(0,1)
    ), data=Minor_Grains2.list , chains=4
    )

#saveRDS(m_mg_reg, "stan_models/m_mg_reg.rds")

Figure 7.10: Bayesian probability estimates for minor grains (Emmer, Einkorn, Foxtail and Broomcorn millets, Sorghum), divided by macroregion and chronology. The prediction estimates are produced with a GLM and weakly informative priors. Credibility intervals are provided.

7.5 Distribution of plants in different site types

I might merge all the fruits (excluding grape and olive) into one graph.

It is possible to examine the distribution of plants across different site types during the four phases under examination. During the Roman period, wheat and barley were most prevalent on rural and religious sites, in the latter being used as part of ritual offerings. Minor grains were widely distributed on religious and urban sites, with rural sites and villas also showing high percentages. Legumes were also used in religious offerings but were found on both urban and rural sites. Grapes were ubiquitous across all Roman sites, with particularly high percentages in urban contexts and being extensively cultivated in large Roman estates. Olives had lower percentages, as they cannot be grown in Northern Italy, where most of the samples were collected. However, they were still widely distributed across all contexts, with peak percentages in urban and religious sites. In the Late Roman period, olive presence significantly decreased in urban sites and villas, but remained relatively unchanged in rural contexts. Similarly to olives, grape decrease in the same contexts, but slightly increase in rural sites. In addition to rural sites, grapes are also diffused in over 66% of the fortified sites (castra) that start to appear in this period. It is noteworthy that fortified sites from the Late Roman period exhibit high percentages of both minor and noble grains, possibly being stocked with the addition of pulses such as faba beans and vetch. In the EMA phase, noble/minor grains and pulses are still present in almost any fortified sites, a situation that is similar to the one in rural sites. In rural sites however minor grains (40%) and pulses (32%) are less ubiquitous in this period, although this figure might suffer from sites that only published cereal remains. Grapes are present in every religious/funerary context, followed by fortified (71%) and rural sites (66%). If olive cultivation seem to decrease in this phase, it might be because most of the sites with archaeobotanical remains are located in Northern Italy. Olives are present in 50% of the religious sites and on 36% of villa sites. Olives are absent from fortified sites in all phases, probably as a reflection of the northern location of these sites. Both rural and urban sites increase the consumption of berries in this period (13-14% ca.), which are also used for funerary offerings (14%).

(a) Noble grains (wheat and barley).

(b) Minor grains.

(c) Legumes.

(d) Grape.

(e) Olive.

(f) Noble grains (wheat and barley).

(g) Minor grains.

Figure 7.11: Ubiquity (%) of plants, divided by site type and chronology.

7.5.1 Bayesian alternative

7.5.1.1 Pre-processing and models

Code

# Eval is set to false because the models have been saved

####################
## NOBLE GRAINS
####################

Noble_Grains <- with(
  Df_Cond_Plants.pa,
  data.frame(
  Chronology = Chronology,
  Macroregion = name_macroreg,
  Site_Type = as.factor(Type),
  Geo = Geo,
  Grain_Count = rowSums(Df_Cond_Plants.pa[,c(11,15)])
  )
)

#Simplify categories
Noble_Grains$Site_Type <- str_replace(Noble_Grains$Site_Type, "Religious, monastery", "Religious")    
Noble_Grains$Site_Type <- str_replace(Noble_Grains$Site_Type, "Castle", "Fortified")    
Noble_Grains$Site_Type <- str_replace(Noble_Grains$Site_Type, "Castrum", "Fortified")  

# Convert to list
Noble_Grain.list <- list(
    NG_Found = Noble_Grains$Grain_Count,
    NG_Tot = 2,
    TC_ID = as.factor((interaction(Noble_Grains$Chronology, 
                                   Noble_Grains$Site_Type)))
   )

Noble_Grain.list$TC_ID <- droplevels(Noble_Grain.list$TC_ID)

# Create a binomial model, where N is 2 because that is the total number of noble 
# grains studied in this thesis
m_ng_type <- ulam(
   alist(
        NG_Found ~ dbinom( 2 , p ),
        logit(p) <- TypeChr[TC_ID],
        TypeChr[TC_ID] ~ dnorm(0,1)
    ), data=Noble_Grain.list , chains=4
    )

####################
## MINOR GRAINS
####################

# Convert to list
Minor_Grains_Type.list <- list(
    MG_Found = Minor_Grains$Grain_Count,
    MG_Tot = 7,
    TC_ID = as.factor((interaction(Minor_Grains$Chronology, 
                                   Minor_Grains$Site_Type)))
   )

Minor_Grains_Type.list$TC_ID <- droplevels(Minor_Grains_Type.list$TC_ID)
  
# Create a binomial model, where N is 2 because that is the total number of noble 
# grains studied in this thesis
m_mg_type <- ulam(
    alist(
        MG_Found ~ dbinom( 7 , p ),
        logit(p) <- TypeChr[TC_ID],
        TypeChr[TC_ID] ~ dnorm(0,1)
    ), data=Minor_Grains_Type.list , chains=4
    )

#saveRDS(m_mg_type, "stan_models/m_mg_type.rds")


####################
## LEGUMES 
####################

Legumes <- with(
  Df_Cond_Plants.pa,
  data.frame(
  Chronology = Chronology,
  Macroregion = name_macroreg,
  Site_Type = as.factor(Type),
  Geo = Geo,
  Seed_Count = rowSums(Df_Cond_Plants.pa[,c(21,27)]) # Excludes Unsp. Pulses
  )
)

#Simplify categories
Legumes$Site_Type <- str_replace(Legumes$Site_Type, "Religious, monastery", "Religious")    
Legumes$Site_Type <- str_replace(Legumes$Site_Type, "Castle", "Fortified")    
Legumes$Site_Type <- str_replace(Legumes$Site_Type, "Castrum", "Fortified")  

# Convert to list
Legumes.list <- list(
    Leg_Found = Legumes$Seed_Count,
    Leg_Tot = 7,
    TC_ID = as.factor((interaction(Legumes$Chronology, 
                                   Legumes$Site_Type)))
   )

Legumes.list$TC_ID <- droplevels(Legumes.list$TC_ID)
  
# Create a binomial model, where N is 2 because that is the total number of noble 
# grains studied in this thesis
m_leg_type <- ulam(
    alist(
        Leg_Found ~ dbinom( 7 , p ),
        logit(p) <- TypeChr[TC_ID],
        TypeChr[TC_ID] ~ dnorm(0,1)
    ), data=Legumes.list , chains=4
    )

#saveRDS(m_leg_type, "stan_models/m_leg_type.rds")

####################
## GRAPE 
####################

Grape <- with(
  Df_Cond_Plants.pa,
  data.frame(
  Chronology = Chronology,
  Macroregion = name_macroreg,
  Site_Type = as.factor(Type),
  Geo = Geo,
  Present = ifelse(Grape>0, 1, 0)
  )
)

#Simplify categories
Grape$Site_Type <- str_replace(Grape$Site_Type, "Religious, monastery", "Religious")    
Grape$Site_Type <- str_replace(Grape$Site_Type, "Castle", "Fortified")    
Grape$Site_Type <- str_replace(Grape$Site_Type, "Castrum", "Fortified")  

# Convert to list
Grape.list <- list(
    Grape_Found = Grape$Present,
    TC_ID = as.factor((interaction(Legumes$Chronology, 
                                   Legumes$Site_Type)))
   )

Grape.list$TC_ID <- droplevels(Grape.list$TC_ID)
  
# Create a binomial model, where N is 2 because that is the total number of noble 
# grains studied in this thesis
m_grape_type <- ulam(
    alist(
        Grape_Found ~ dbinom( 1 , p ), # Basically a bernoulli dist
        logit(p) <- TypeChr[TC_ID],
        TypeChr[TC_ID] ~ dnorm(0,2)
    ), data=Grape.list , chains=4
    )

#saveRDS(m_grape_type, "stan_models/m_grape_type.rds")

####################
## OLIVES 
####################

Olives <- with(
  Df_Cond_Plants.pa,
  data.frame(
  Chronology = Chronology,
  Macroregion = name_macroreg,
  Site_Type = as.factor(Type),
  Geo = Geo,
  Present = ifelse(Olive>0, 1, 0)
  )
)

#Simplify categories
Olives$Site_Type <- str_replace(Olives$Site_Type, "Religious, monastery", "Religious")
Olives$Site_Type <- str_replace(Olives$Site_Type, "Castle", "Fortified")    
Olives$Site_Type <- str_replace(Olives$Site_Type, "Castrum", "Fortified")  

# Convert to list
Olive.list <- list(
    Olive_Found = Olives$Present,
    TC_ID = as.factor((interaction(Olives$Chronology, 
                                   Olives$Site_Type)))
   )

Olive.list$TC_ID <- droplevels(Olive.list$TC_ID)
  
# Create a binomial model, where N is 2 because that is the total number of noble 
# grains studied in this thesis
m_olives_type <- ulam(
    alist(
        Olive_Found ~ dbinom( 1 , p ), # Basically a bernoulli dist
        logit(p) <- TypeChr[TC_ID],
        TypeChr[TC_ID] ~ dnorm(0,2)
    ), data=Olive.list , chains=4
    )

#saveRDS(m_olives_type, "stan_models/m_olives_type.rds")

####################
## OLIVES 
####################

Nuts <- with(
  Df_Cond_Plants.pa,
  data.frame(
  Chronology = Chronology,
  Macroregion = name_macroreg,
  Site_Type = as.factor(Type),
  Geo = Geo,
  Seed_Count = rowSums(Df_Cond_Plants.pa[,c(29,31)])
  )
)

#Simplify categories
Nuts$Site_Type <- str_replace(Nuts$Site_Type, "Religious, monastery", "Religious")    
Nuts$Site_Type <- str_replace(Nuts$Site_Type, "Castle", "Fortified")    
Nuts$Site_Type <- str_replace(Nuts$Site_Type, "Castrum", "Fortified")  

# Convert to list
Nuts.list <- list(
    Nuts_Found = Nuts$Seed_Count,
    Nuts_Tot = 3,
    TC_ID = as.factor((interaction(Nuts$Chronology, 
                                   Nuts$Site_Type)))
   )

Nuts.list$TC_ID <- droplevels(Nuts.list$TC_ID)
  
# Create a binomial model, where N is 3 because that is the total number of noble 
# grains studied in this thesis
m_nuts_type <- ulam(
    alist(
        Nuts_Found ~ dbinom( 3 , p ),
        logit(p) <- TypeChr[TC_ID],
        TypeChr[TC_ID] ~ dnorm(0,1)
    ), data=Nuts.list , chains=4
    )


#saveRDS(m_nuts_type, "stan_models/m_nuts_type.rds")

7.5.1.2 Results

(a) Noble grains (wheat and barley).

(b) Minor grains.

(c) Legumes.

(d) Grape.

(e) Olives. The model is predicting (conservative) values for fortified sites, in which no olive was found.

(f) Nuts.

(g) Domestic fruits.

Figure 7.12: Bayesian probability estimates for main plant taxa, divided by site type and chronology. The prediction estimates are produced with a GLM and weakly informative priors. Credibility intervals are provided.

7.5.2 Richness and diversity in the sites

(a) Roman plant richness

(b) Roman plant diversity.

(c) Late Roman plant richness

(d) Late Roman plant diversity.

(e) EMA plant richness

(f) EMA plant diversity.

Figure 7.13: Violin plots of plant richness and diversity in different site types. The grey dots (jitters) indicate the value for the single site, while the white boxplot shows the median and the quartile values.

7.5.2.1 Richness in urban contexts

The richness of plants in urban contexts has been calculated for every phase. As not every sample reported cereals, fruits/nuts and legumes, a correction has also been applied to exclude contexts that did not report the three types of plant taxa under examination. This correction has however strong disadvantages as it removes sites where the absence of certain plant categories may have been deliberate (e.g., storages). However, if the figures vary considerably (especially in the case of the 11^th century CE which only consists of 2 samples after the correction), the general trend is similar after the correction. Plant richness in urban contexts slightly decreases after the Roman age, but early medieval cities are richer in plants when compared to the Roman and Late Roman phases. In the 11th century, plant richness increases further. These values will be discussed in the results chapter.

Table 7.2: Plant richness in urban contexts before and after (*) removing incomplete samples.

Chronology	Richness	Samples	Richness*	Samples*
R	6.5	44	8.45	22
LR	5.9	24	7.2	10
EMA	7.2	25	9.6	8
Ma	9	8	17.5	2

7.5.3 Correspondence analysis tests

7.6 Soil types vs Crops

see Rippon et al 2015, p. 81 as a comparison for results and Lodwick chapter in Hamerow and McKerracher 2022, p. 153

8 To do:

8.1 Wine/Olive oil production in Italy

Questions: Is there a decrease in the production of cash crops in Late Roman villas? In particular in central Italian villas.

Available data:

• Seed remains

• Wine presses in Roman Italy (1st to 6th c. ca) -> Dodd

• Wine/Olive presses, dated but no info on wine/olive or size: http://oxrep.classics.ox.ac.uk/databases/olive_oil_and_wine_presses_database/

• Dodd: “a notable increase in size from around the first century BCE.” No quantitative data provided.

8.2 Fruit consumption

• The domesticated and the wild

• What fruits arrive on the Roman table?

• Can we see a decreased consumption of figs during the LALIA in the north?

This kind of data manipulation might be misleading as it is ignoring the fact that some centuries are only represented by - let’s say - two sites. Is this a good idea?

Figure 8.1: The values for Central Italy and for Late Roman Southern Italy are based on very few contexts.

Anderson, M.J., 2017. Permutational Multivariate Analysis of Variance (PERMANOVA). In: Wiley StatsRef: Statistics Reference Online. John Wiley & Sons, Ltd, pp. 1–15.

Bowes, K.D., University of Pennsylvania (Eds.), 2020. The Roman Peasant Project 2009-2014: Excavating the Roman rural poor, First edition. ed, University museum monograph. University of Pennsylvania, Museum of Archaeology and Anthropology, Philadelphia.

Kosub, S., 2019. A note on the triangle inequality for the Jaccard distance. Pattern Recognition Letters 120, 36–38.

Oksanen, J., Blanchet, F.G., Friendly, M., Kindt, R., Legendre, P., McGlinn, D., Minchin, P.R., O’Hara, R.B., Simpson, G.L., Solymos, P., Stevens, M.H.H., Szoecs, E., Wagner, H., 2020. Vegan: Community ecology package.

---

1. The Late Roman values for Southern Italy are only based on 5 samples (3 of which are from the same site, Salapia) so the values are not very trustworthy.

2. Source: Wikipedia. Change the source later